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Abstract 

Annealed importance sampling is a means to assign equilibrium 
weights to a nonequilibrium sample that was generated by a simu- 
lated annealing protocol pp. The weights may then be used to calculate 
equilibrium averages, and also serve as an "adiabatic signature" of the 
chosen cooling schedule. In this paper we demonstrate the method on 
the 50-atom dileucine peptide, showing that equilibrium distributions 
are attained for manageable cooling schedules. For this system, as 
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naively implemented here, the method is modestly more efficient than 
constant temperature simulation. However, the method is worth con- 
sidering whenever any simulated heating or cooling is performed (as is 
often done at the beginning of a simulation project, or during an NMR 
structure calculation), as it is simple to implement and requires mini- 
mal additional CPU expense. Furthermore, the naive implementation 
presented here can be improved. 

1 Introduction 

Simulated annealing (SA) is used in a wide variety of biomolecular 
calculations. Crystallographic refinement protocols [2] and standard 
NMR structure calculations [3j [U [5] both rely on SA to optimize a 
"target function," constructed so that the global minimum of the tar- 
get function corresponds to the native structure. Molecular dynamics 
calculations often begin by cooling a configuration from a high tem- 
perature ensemble to a lower temperature, at which the simulation is 
to be performed. 

In this paper, we consider a different use for SA calculations. Since 
a set of structures that is generated by a series of SA trajectories is a 
nonequilibrium sample, they may not be used to calculate equilibrium 
averages. However, Neal demonstrated a simple procedure, called "an- 
nealed importance sampling" (AIS) that allows the nonequilibrium 
sample to be reweighted into an equilibrium one[lJ. AIS is closely 
connected with the Jarzynski relation[6j. To our knowledge, the al- 
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gorithm has only appeared once in the chemical physics literature [7], 
where it was used (along with sophisticated Monte Carlo techniques) 
to sample a one-dimensional potential. Here, we demonstrate an ap- 
plication of the AIS algorithm to generate an equilibrium sample of 
an implicitly solvated peptide, and discuss other uses for AIS which 
may of interest to the molecular simulation community. 

The basic idea which underlies SA is also the motivation for other 
temperature based sampling methods, notably J-walking[8j, simulated 
tempering^ [10] and replica exchange/parallel tempering[TT| [T2]. By 
coupling a simulation to a high temperature reservoir, it is hoped that 
the low temperature simulation may explore the configuration space 
more thoroughly. This is achieved by thermally activated crossing 
of energetic barriers, which are large compared to the thermal en- 
ergy scale of the lower temperature simulation, but are crossed more 
frequently at higher temperature. Simulated and parallel tempering 
differ in the way that the different temperature simulations are cou- 
pled. Simulated tempering heats and then cools the system, in a way 
that maintains an equilibrium distribution. Parallel tempering couples 
simulations run in parallel at different temperatures by occasionally 
swapping configuartions between temperatures, again in such a way 
that canonical sampling is maintained. 

AIS offers yet another approach to utilizing a high temperature 
ensemble for equilibrium sampling at a lower temperature. A sam- 
ple of a high temperature ensemble is annealed to a lower temper- 
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ature, by alternating constant temperature simulation with steps in 
which the tempertaure is jumped to a lower value. Each annealed 
structure is assigned a weight, which depends on the trajectory that 
was traced during the annealing process. Equilibrium averages over 
the lower temperature ensemble may then be calculated by a simple 
weighted average. Furthermore, the distribution of trajectory weights 
contains useful information about the statistics of the annealed sam- 
ple. Roughly, a schedule which quenches high temperature structures 
very rapidly to low temperature will result in a sample dominated by 
a few high weight structures, resulting in poor statistics. This con- 
nection between the distribution of weights and the extent to which 
the schedule is not adiabatic ought to be of interest to anyone who 
uses SA protocols — whether for equilibrium sampling or for structure 
calculation. 

We have used the AIS method to generate 298 K equilibrium en- 
sembles of the dileucine peptide, by annealing structures from a 500 
K distribution with several different cooling schedules. For the most 
efficient schedule used, we found a modest gain (about a factor of 3) 
over constant temperature simulation. This result is consistent with 
earlier observations on the expected efficiency of temperature-based 
sampling methods [13]. 
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2 Theory 



Consider a standard simulated annealing (SA) trajectory, in which a 
protein is slowly cooled from a conformation x at a (high) temperature 
Tm- The cooling is achieved by alternating constant temperature 
dynamics with "temperature jumps," during which the temperature 
is lowered instantaneously. Usually, the system is cooled to a low 
temperature, since the aim of standard SA calculations is to find the 
global minimum on the energy landscape. But we can imagine instead 
ending the run at To = 300 K — in fact, we can think of many such 
runs, all ending at 300 K. We then have an ensemble of conformations, 
though clearly not distributed canonically at To. We would like to 
know if there is a way to reweight this distribution, so that it can be 
used to compute equilibrium averages at To. The affirmative answer 
is provided by the annealed importance sampling (AIS) method. 

To make the discussion more concrete, consider many independent 
annealing trajectories Xj(t) which at time tu-i have just been cooled 
from inverse temperature (5m to (3m -i- As usual, each temperature 
defines a distribution of conformations: 7Tj(x) oc exp[— /3jC7(x)]. Imme- 
diately after £m-i, before the system is allowed to relax to 7Tm-i(x), 
we can compute the equilibrium average of an arbitrary quantity A 
over 7Tm-i(x) by using the weight u>(x) = 7Tm_i(x)/7Tm(x): 

(A) M -lZ M -l = / dx^(x)7T M -l(x) = / dx^(x)7T M (x)w(x), (1) 
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where (A)i denotes an average over 7Tj, and Zi = j dx7Tj(x). In other 
words, we may reweight the distribution 7rj\/(x) to calculate averages 
over 7Tm_i(x), by multiplying by the ratio of Boltzmann factors. 

Generalizing the argument to M temperature steps is straightforward PQ, 
by forming the product of weights for successive cooling steps: 

I 1+ \\ TT 7r *— 1 ( X J 1 )) / n 

^- = ^(to))=n Ji(x - (ti , l)) ■ p) 

Equation [2] gives the weight for trajectory j, cooled at successive times 
tM-i, tM-2,--- through inverse temperatures (3m, Pm-i,--- to reach 
conformation x,-(io)- At each temperature, reweighting ensures that 
averages may be calculated for the appropriate canonical distribution, 
even though the system has not yet relaxed. 

The AIS idea is easily turned into an algorithm for producing a 

canonical distribution from serially generated annealing trajectories: 

(i) Generate a sample of the distribution ttm(x), by a 

sufficiently long simulation at Tm- 

(ii) Pull a conformation from 7Tm( x ) at random and anneal 
down to /3o, yielding conformation xi(to). Keep track of the 
weight io(xi(io)) for this trajectory by Eq. 

(iii) Repeat steps (hi) and (iv) N times, yielding congiura- 
tions Xj and weights w(xj) = Wj for j = 1, 1, N. 

Equilibrium averages at temperature Tq are then calculated by a 
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weighted average: 

(A) = ^ 3 3 (3) 

The cooling schedule is defined by the number and spacing of the 
temperature steps, as well as the duration of the constant temperature 
simulation at each step. As available resources necessarily limit the 
CPU time spent on each annealing trajectory, careful consideration of 
the schedule is in order. Clearly, a schedule in which high temperature 
configurations are quenched in one step to low temperature amounts 
to a single-step reweighting procedure[14j. We may expect that such a 
schedule would be quite ineffective for large temperature jumps, since 
very few configurations in the high temperature distribution have ap- 
preciable weight in the low temperature distribution. By introducing 
intermediate steps, the system is allowed to relax locally, bridging the 
high and low temperature distributions in a way that echoes replica 
exchange protocols |ll[|12j . simulated tempering [£l [TO], and the multi- 
ple histogram method [TO]. However, the "top-down" structure of the 
algorithm most closely resembles J- walking |8l 116] . 

3 Results 

The dileucine peptide (ACE-[Leu]2-NME) is good choice for the vali- 
dation of new algorithms, as it is small enough (50 atoms, including 
nonpolar hydrogens) that exhaustive sampling by standard simulation 
methods is possible, yet more akin to protein systems than a one- or 
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two-dimensional "toy" model. 

The high temperature ensemble was generated by 300 nsec of 
Langevin dynamics at Tm = 500 K, as implemented in Tinker v. 
4.2.2 [17J , with a timestep of 1.0 fsec, and a friction constant of 91 
psec -1 , and solvation was treated by the GB/SA method [IS]. Frames 
were written every psec, resulting in a sample of 3 x 10 4 frames in the 
high temperature sample. 

The 500 K sample was annealed down to 298 K using 4 different 
schedules, consisting of a total of 3, 5, 9, and 17 temperatures, includ- 
ing the endpoints. In each case, the temperatures were distributed 
geometrically. Following each temperature jump, the velocities were 
reinitialized by sampling randomly from the Maxwell-Boltzmann dis- 
tribution, and then allowed to relax at constant temperature for a 
time tn = 0.5 psec (except where noted) with the protocol described 
above. A total of A = 1.6 x 10 4 annealing trajectories were generated 
for each schedule. The control of the integration routine to effect the 
annealing, as well as the calculation of the trajectory weights, were 
implemented in a Perl script. 

Figure [I] shows that the 298 K distribution of energy is recovered 
by the AIS procedure. It is noteworthy that the 500 K distribution 
(corresponding to the high T sample) overlaps very little with the 298 
K distribution, and yet the 298 K distribution is reproduced well for 
the two slowest schedules. Equally interesting is how poorly the algo- 
rithm performs when the structures are cooled too rapidly, especially 
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on the low E side of the distribution, where there is no overlap with 
the high T distribution. We conclude that the schedules with 3 or 5 
T-steps quench the structures too rapidly, resulting in many of the tra- 
jectories becoming "stuck" in high-energy states that are metastable 
at 298 K. 

This last observation may be quantified by asking, "How many 
of the annealed structures contribute appreciable weight to averages 
calculated with Eq. [3p" To address this question, for each schedule we 
estimated the number of configurations n which contribute appreciable 
weight to the averages: 



n = ^lH = /tf, (4) 

^ max 



where w max is the largest weight observed (see Table[T]). If this number 
is near 1, then a small number of trajectories dominate the average — 
see Eq.[3] — and poor results should be expected. The effective fraction 
of the annealing trajectories which generate "useful" or "successful" 
structures is denoteed by /. 

A more complete picture is provided by the full distribution of 
the (logarithm of) trajectory weights (Fig. [2]). For each schedule, the 
weights which contribute the most to the T = 298 K sample are to 
the right, at large values of w. The trend is clear — as slower cooling 
is effected, the distribution narrows and shifts to the right. It has 
been shown that the accuracy of averages computed from this type of 
protocol is roughly related to the variance of the (adjusted) weights[T]. 



9 



(The adjusted weight is the weight divided by the average weight.) 
This "rule of thumb" is borne out by the data in Fig. [2] and Table 
[I] — as the cooling slows down the distribution of weights narrows, and 
the number of trajectories contributing to the equilibrium averages 
increases. This type of analysis may serve as a means of distinguishing 
between annealing schedules to decide on a cooling schedule which is 
slow enough to yield reasonable estimates of equilibrium averages. It 
is also essential for optimizing an AIS protocol for sampling efficiency, 
as discussed in the next few paragraphs. 

How much better than standard simulation (if at all) is equilib- 
rium sampling by AIS? In order to make a direct comparison between 
AIS and constant temperature simulation, we need to compare the 
CPU time invested per statistically independent configuration in each 
protocol. For the constant temperature simulation, this time may be 
estimated in several ways[19j [20], and is essentially the time needed 
for the simulation to "forget" where it has been. Following the con- 
vention for correlation times, we call this time r, = t(Tj), where i 
labels the temperature: M for the high T distribution, and for the 
low T distribution. For the system studied here, tm = 0.8 nsec and 
ro = 3.0 nsec, as estimated from timseries of the a — ► (3 backbone 
dihedral transition [13]. 

The total cost to generate a structure in an AIS simulation is the 
sum of the costs of generating a structure in the high T distribution 
plus that for the annealing phase. Of course, not every annealing 
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trajectory contributes to thermodynamic averages (Eq. [3|). What then 
is the total cost t cos t of a "successful" annealed structure? The first 
part is from high temperature sampling — i.e., tm- The second part is 
the cost of all the annealing trajectories, divided by the number which 
contribute to equilibrium averages. The time ianneal is the time spent 
annealing each structure: 



Recall that tn is the duration of the constant temperature relaxation 
steps, and there is no relaxation phase at the highest and lowest tem- 
peratures. 

The total cost t cos t is then the sum of tm and tanneai: 



The efficiency of an AIS protocol may then be computed by taking 
the ratio R = tq/ t cost (see Table [1]) , which gives the factor by which 
an AIS protocol is more or less efficient than constant temperature 
simulation. The data in Table [1] show that the best schedule used 
here offer a modest speedup over constant temperature simulation, of a 
factor of about 3. These findings are in agreement with an analysis we 
have published of another temperature-based sampling protocol[13j. 
We note that an optimized AIS protocol would require tuning N based 
on (perhaps preliminary) estimates of /. 



^anneal = tn(M - 2) 



(5) 



tc — T M + ^anneal/ /■ 



(6) 



11 



It is instructive to compare the AIS results to simple reweighting — 
i.e., AIS with no intermediate temperature steps or relaxation. In 
this case, no computer time is spent annealing, and the efficiency 
gain is simply tq/tm = 3.75. The fraction / is of course reduced 
compared to any AIS protocol — when reweighting our 500 K dileucine 
trajectory to 298 K distribution, / = 1.3 x 10~ 4 — but this has no 
impact on the efficiency, provided a sufficient number of snapshots are 
available for reweighting. However, it is clear that / will be greatly 
reduced in systems which undergo a folding transition upon lowering 
the temperature. This is simply a reflection of the fact that there is 
negligible overlap between the folded and unfolded distributions. In 
such useful reweighting protocol would require the generation 

of astronomical numbers of structures in the Tm distribution, and 
annealing is advised. 

4 Conclusion 

We have demonstrated the application of Neal's annealed importance 
sampling (AIS) algorithm for equilibrium sampling of the dileucine 
peptide. AIS allows the calculation of equilibrium averages from a 
nonequilibrium sample of strutures that results from a simulated an- 
nealing protocol. To our knowledge, AIS has not previously been 
applied to a molecular system. While the method, as naively im- 
plemented here, represents only a modest improvement over constant 
temperature simulation, it is interesting for several reasons beyond 
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equilibrium sampling. 

First, in applications where simulated annealing is already in widespread 
use (most notably, NMR structure calculations [3j HI [5]), the path 
weights may be used to calculate (perhaps noisy) equilibrium aver- 
ages, and perhaps ultimately Boltzmann-distributed ensembles. The 
path weights also contain information that can be used to discriminate 
between different schedules, which may provide a way to optimize the 
schedule, based on the analysis of t a nn> the cost of annealing to "good" 
structures. 

Second, it may be possible to improve considerably on the efficiency 
of the method by implementing a more sophisticated version, which 
uses a resampling procedure to prune the low weight paths at each 
cooling step. (For a detailed discussion of resampling methods, see 
the book by Liu[2T].) In this approach, we first cool some number N 
of structures from the high temperature (Tm) ensemble, yielding N 
weighted structures at Tm-i- We then resample ./V times from this 
Tm-i ensemble, according to the cumulative distribution function of 
the weights, pruning the low weight paths without biasing the sample. 
This type of approach was recently applied successfully to sampling 
near native protein configurations of a discretized and coarse-grained 
model [22]. Nevertheless, we emphasize that the ultimate efficiency of 
any AIS protocol limited by the intrinsic sampling rate of the highest 
temperature, which may be modest; see Ref. ??. 

Finally, the AIS procedure could be naturally combined with "an- 
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nealing" in the parameters of the Hamiltonian. Such a hybrid of AIS 
and Hamiltonian switching might be used, for example, to transform 
an NMR target function into a molecular mechanics potential function, 
over the course of a structure calculation. The result of such a cal- 
culation would be an equilibrium ensemble of structures, distributed 
according to the molecular mechanics potential. Such ensembles would 
find wide application, for instance in docking or homology modeling. 
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Table 1: Comparison of the efficiency of AIS between several cooling sched- 



ules, n is given by Eq. 0], t cost is given by Eq. [6] The efficiency gain is the total 
simulation time invested in each successful annealed structure t cost divided 
by the time needed to generate an indepenendent structure by constant tem- 
perature simulation. The t indicates schedules for which data are presented 
in Figs. [1] and [2j 
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Figure Legends 

Figure [H 

Distribution of energies, from standard, constant temperature simu- 
lation and AIS. The dashed line is the T = 500 K distribution that 
was used for the high T ensemble. The other data compare a 300 
nsec, T = 298 K constant temperature simulation to 298 K ensembles 
generated by the AIS algorithm with different cooling schedules. The 
schedules are discussed in Table [TJ 

Figure [21 

Distribution of the logarithm of trajectory weights for the four cooling 
schedules used in Fig. [1] and discussed in Table [TJ 
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